Figures: Conditional Logistic Results

1. Load results & prepare data table for plots

load(res_fname, verbose = T)

## Loading objects:
##   res_alg
##   res_arh
##   res_bmi
##   res_crd
##   res_crd_wbmi
##   res_dir
##   res_fa
##   res_mets_bd
##   res_mets_bd_bmi
##   res_mets_ics
##   res_mets_ics_bmi
##   res_mets_ocs
##   res_mets_ocs_bmi
##   res_normal
##   res_red_mlk
##   mets_info
##   mets_pheno
##   mets_pheno_asth
##   mets_pheno_trim
##   cer_list
##   ceroth_list
##   sm_list
##   sph_list
##   sphoth_list

## Model results to show

res_long <- rbind(res_crd[, model := "crude"], 
                  res_bmi[, model := "BMI adjusted (n-pairs=293)"], 
                  res_crd_wbmi[, model := "in subject w BMI (not BMI adjusted)"], 
                  res_normal[, model := "exclude red/milky samples (n-pairs=545)"])
res_long[, dir_eff := case_when(or > 1 ~ "positive", 
                                or <= 1 ~ "negative")]
res_long[, dir_eff := factor(dir_eff, levels = c("negative", "positive"))]
res_long[, model := factor(model, levels = c("crude", "BMI adjusted (n-pairs=293)", "in subject w BMI (not BMI adjusted)", 
                                             "exclude red/milky samples (n-pairs=545)"))]
res_long[, metabolite := factor(metabolite, levels = c(sphoth_list, sm_list, cer_list, ceroth_list))]
res_long[, rev_metabolite := factor(metabolite, levels = rev(levels(metabolite)))]


## Models results with minimal change from crude

res_unused <- rbind(res_crd[, model := "crude"], 
                    res_arh[, model := "allergic rhinitis adjusted"], 
                    res_fa[, model := "food allergy #Dx adjusted"], 
                    res_alg[, model := "unspecified allergy #Dx adjusted"])
res_unused[, dir_eff := case_when(or > 1 ~ "positive", 
                                  or <= 1 ~ "negative")]
res_unused[, dir_eff := factor(dir_eff, levels = c("negative", "positive"))]
res_unused[, model := factor(model, levels = c("crude", "allergic rhinitis adjusted", "food allergy #Dx adjusted", 
                                               "unspecified allergy #Dx adjusted"))]
res_unused[, metabolite := factor(metabolite, levels = c(sphoth_list, sm_list, cer_list, ceroth_list))]
res_unused[, rev_metabolite := factor(metabolite, levels = rev(levels(metabolite)))]

1. Volcano plots

1.1. Crude model results

res_crd[fdr_bh < alpha_thld, max(pval)]

## [1] 0.01291851

res_crd[fdr_bh >= alpha_thld, min(pval)]

## [1] 0.01539316

fdr_thld <- 0.014

bonf_thld <- 0.05/mets_info[, length(metabolite)]

ggplot(res_crd, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "black", linetype = 2) + 
        geom_hline(yintercept = -log10(fdr_thld), color = "#00798c", linetype = 2) + 
        geom_hline(yintercept = -log10(bonf_thld), color = "#edae49", linetype = 2) + 
        labs(title = "Crude conditional logistic model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = 0, color = "black", size = 4) + 
        annotate("text", x = 0, y = -log10(fdr_thld), label = "FDR_BH == 0.05", parse = T, 
                 vjust = 1.5, hjust = 0, color = "#00798c", size = 4) + 
        annotate("text", x = 0, y = -log10(bonf_thld), label = "Bonferroni threshold", parse = F, 
                 vjust = 1.5, hjust = 0, color = "#edae49", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < fdr_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_crd.jpg"), width = 12, height = 7)

ggplot(res_crd_wbmi, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Crude conditional logistic model results in subjects with BMI (not BMI-adjusted)", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = 0, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_crd_subj_wbmi.jpg"), width = 12, height = 7)

1.2. Restricted to normal blood samples

res_normal[fdr_bh < alpha_thld, max(pval)]

## [1] 0.009423959

res_normal[fdr_bh >= alpha_thld, min(pval)]

## [1] 0.01132697

fdr_thld <- 0.01

ggplot(res_normal, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "black", linetype = 2) + 
        geom_hline(yintercept = -log10(fdr_thld), color = "#00798c", linetype = 2) + 
        geom_hline(yintercept = -log10(bonf_thld), color = "#edae49", linetype = 2) + 
        labs(title = "Crude conditional logistic model results - exclude red/milky samples", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = 0, color = "black", size = 4) + 
        annotate("text", x = 0, y = -log10(fdr_thld), label = "FDR_BH == 0.05", parse = T, 
                 vjust = 1.5, hjust = 0, color = "#00798c", size = 4) + 
        annotate("text", x = 0, y = -log10(bonf_thld), label = "Bonferroni threshold", parse = F, 
                 vjust = 1.5, hjust = 0, color = "#edae49", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < fdr_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_normal.jpg"), width = 12, height = 7)

1.3. BMI adjusted results

ggplot(res_bmi, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "BMI-adjusted conditional logistic model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_bmi.jpg"), width = 12, height = 7)

1.4. Association with ICS/OCS

## ICS

ggplot(mets_pheno_trim[Group == "Asthmatic", ]) + 
        geom_histogram(aes(inhaled_steriod_prescriptions), bins = 30) + 
        labs(title = "Number of ICS prescriptions in asthmatics") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "ics_asthmatics.jpg"), width = 12, height = 7)

ggplot(res_mets_ics, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of ICS prescriptions - linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_ics.jpg"), width = 12, height = 7)

ggplot(res_mets_ics_bmi, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of ICS prescriptions - BMI-adjusted linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld*2, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_ics_bmi.jpg"), width = 12, height = 7)


## OCS

ggplot(mets_pheno_trim[Group == "Asthmatic", ]) + 
        geom_histogram(aes(oral_steroid_prescriptions), bins = 30) + 
        labs(title = "Number of OCS prescriptions in asthmatics") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "ocs_asthmatics.jpg"), width = 12, height = 7)

ggplot(res_mets_ocs, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of OCS prescriptions - linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_ocs.jpg"), width = 12, height = 7)

ggplot(res_mets_ocs_bmi, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of OCS prescriptions - BMI-adjusted linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld*2, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_ocs_bmi.jpg"), width = 12, height = 7)


## Bronchodilator

ggplot(mets_pheno_trim[Group == "Asthmatic", ]) + 
        geom_histogram(aes(bronchodialator_prescriptions), bins = 30) + 
        labs(title = "Number of bronchodilator prescriptions in asthmatics") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "bd_asthmatics.jpg"), width = 12, height = 7)

ggplot(res_mets_bd, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of bronchodilator prescriptions - linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_bd.jpg"), width = 12, height = 7)

ggplot(res_mets_bd_bmi, aes(x = beta, y = -log10(pval))) + 
        geom_point(alpha = 0.4, size = 2.5) + 
        geom_hline(yintercept = -log10(alpha_thld), color = "#00798c", linetype = 2) + 
        labs(title = "Metabolites & number of bronchodilator prescriptions - BMI-adjusted linear model results", 
             x = expression('Estimated '*beta), 
             y = expression('-log'[10]*'('*italic(P)*'-value)')) + 
        annotate("text", x = 0, y = -log10(alpha_thld), label = "P_value == 0.05", parse = T, 
                 vjust = 1.5, hjust = -1.8, color = "#00798c", size = 4) + 
        geom_label_repel(aes(label = ifelse(pval < alpha_thld, metabolite, "")),
                  box.padding   = 0.25, 
                  point.padding = 0.5,
                  segment.color = 'grey50') + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12))

ggsave(here(fig_dir, "volc_mets_bd_bmi.jpg"), width = 12, height = 7)

2. OR & CI plots for metabolites

2.1. Model results to show

2.1.1. Sphingolipids other

ggplot(res_long[metabolite %in% sphoth_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_long[metabolite %in% sphoth_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "sphoth_or95ci.jpg"), width = 14, height = 7)

2.1.2. SM

ggplot(res_long[metabolite %in% sm_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_long[metabolite %in% sm_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "sm_or95ci.jpg"), width = 14, height = 7)

2.1.3. Cermides

ggplot(res_long[metabolite %in% cer_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_long[metabolite %in% cer_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "cer_or95ci.jpg"), width = 14, height = 7)

ggplot(res_long[metabolite %in% ceroth_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_long[metabolite %in% ceroth_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "ceroth_or95ci.jpg"), width = 14, height = 7)

2.2. Model results with minimal change from crude

2.2.1. Sphingolipids other

ggplot(res_unused[metabolite %in% sphoth_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_unused[metabolite %in% sphoth_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "sphoth_or95ci_min_change_crd.jpg"), width = 14, height = 7)

2.2.2. SM

ggplot(res_unused[metabolite %in% sm_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_unused[metabolite %in% sm_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "sm_or95ci_min_change_crd.jpg"), width = 14, height = 7)

2.2.3. Cermides

ggplot(res_unused[metabolite %in% cer_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_unused[metabolite %in% cer_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "cer_or95ci_min_change_crd.jpg"), width = 14, height = 7)

ggplot(res_unused[metabolite %in% ceroth_list, ], aes(x = rev_metabolite, y = or)) + 
        facet_grid(. ~ model) + 
        geom_errorbar(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                      width = 0.25) + 
        geom_pointrange(aes(ymin = lower95, ymax = upper95, color = dir_eff), 
                        size = 0.5) + 
        scale_colour_hp_d(option = "Ravenclaw", name = "direction of effect") + 
        geom_hline(yintercept = 1, linetype = 2) + 
        labs(title = "Metabolites & asthma status", 
             x = "Metabolite", y = "Odds Ratio (95% Confidence Interval)") + 
        theme_minimal() + 
        theme(title = element_text(size = 16), 
              strip.text = element_text(size = 12), 
              axis.title = element_text(size = 12), 
              axis.text = element_text(size = 12), 
              legend.position = "bottom",
              legend.title = element_text(size = 12), 
              legend.text = element_text(size = 12)) + 
        scale_y_continuous(trans = log_trans(), 
                           breaks = axisTicks(log(range(res_unused[metabolite %in% ceroth_list, or])), 
                                              log = TRUE)) + 
        coord_flip() 

ggsave(here(fig_dir, "ceroth_or95ci_min_change_crd.jpg"), width = 14, height = 7)

4. Session info

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS release 6.10 (Final)
## 
## Matrix products: default
## BLAS:   /app/R-3.6.0@i86-rhel6.0/lib64/R/lib/libRblas.so
## LAPACK: /app/R-3.6.0@i86-rhel6.0/lib64/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.1.0      harrypotter_2.1.0 ggrepel_0.8.1     here_0.1         
##  [5] data.table_1.12.8 forcats_0.4.0     stringr_1.4.0     dplyr_0.8.4      
##  [9] purrr_0.3.3       readr_1.3.1       tidyr_1.0.2       tibble_2.1.3     
## [13] ggplot2_3.2.1     tidyverse_1.3.0  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.0.0 xfun_0.12        reshape2_1.4.3   haven_2.2.0     
##  [5] lattice_0.20-38  colorspace_1.4-1 vctrs_0.2.2      generics_0.0.2  
##  [9] htmltools_0.4.0  yaml_2.2.1       rlang_0.4.4      pillar_1.4.3    
## [13] withr_2.1.2      glue_1.3.1       DBI_1.1.0        dbplyr_1.4.2    
## [17] modelr_0.1.5     readxl_1.3.1     plyr_1.8.5       lifecycle_0.1.0 
## [21] munsell_0.5.0    gtable_0.3.0     cellranger_1.1.0 rvest_0.3.5     
## [25] evaluate_0.14    labeling_0.3     knitr_1.28       fansi_0.4.1     
## [29] broom_0.5.4      Rcpp_1.0.3       backports_1.1.5  jsonlite_1.6.1  
## [33] farver_2.0.3     fs_1.3.0         gridExtra_2.3    hms_0.5.3       
## [37] digest_0.6.23    stringi_1.4.5    rprojroot_1.3-2  grid_3.6.0      
## [41] cli_2.0.1        tools_3.6.0      magrittr_1.5     lazyeval_0.2.2  
## [45] crayon_1.3.4     pkgconfig_2.0.3  xml2_1.2.2       reprex_0.3.0    
## [49] lubridate_1.7.4  assertthat_0.2.1 rmarkdown_2.1    httr_1.4.1      
## [53] rstudioapi_0.11  R6_2.4.1         nlme_3.1-144     compiler_3.6.0